import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import plotly.express as px
import seaborn as sns
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
from scipy.stats import chi2_contingency
from statsmodels.stats.contingency_tables import mcnemar
from cmh import CMH
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, median_absolute_error, r2_score
from sklearn.preprocessing import PolynomialFeatures
from outliers import smirnov_grubbs as grubbs
from neulab.OutlierDetection import DixonTest
from scipy.stats import kstest, shapiro, anderson, cramervonmises, norm
from statsmodels.stats.diagnostic import lilliefors
from numpy.random import randn
from scipy.stats import ttest_ind
import statsmodels.stats.power as smp

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import warnings
warnings.filterwarnings('ignore')

Загрузка данных

Данные: Heart Attack dataset - 13 медицинских признаков и колонка о наличии сердечного приступа

Источник: https://www.kaggle.com/datasets/pritsheta/heart-attack

df = pd.read_csv('Heart_Attack.csv')
df.head()

	age	sex	cp	trestbps	chol	fbs	restecg	thalach	exang	oldpeak	slope	ca	thal	target
0	63	1	3	145	233	1	0	150	0	2.3	0	0	1	1
1	37	1	2	130	250	0	1	187	0	3.5	0	0	2	1
2	41	0	1	130	204	0	0	172	0	1.4	2	0	2	1
3	56	1	1	120	236	0	1	178	0	0.8	2	0	2	1
4	57	0	0	120	354	0	1	163	1	0.6	2	0	2	1

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 303 entries, 0 to 302
Data columns (total 14 columns):
 #   Column    Non-Null Count  Dtype  
---  ------    --------------  -----  
 0   age       303 non-null    int64  
 1   sex       303 non-null    int64  
 2   cp        303 non-null    int64  
 3   trestbps  303 non-null    int64  
 4   chol      303 non-null    int64  
 5   fbs       303 non-null    int64  
 6   restecg   303 non-null    int64  
 7   thalach   303 non-null    int64  
 8   exang     303 non-null    int64  
 9   oldpeak   303 non-null    float64
 10  slope     303 non-null    int64  
 11  ca        303 non-null    int64  
 12  thal      303 non-null    int64  
 13  target    303 non-null    int64  
dtypes: float64(1), int64(13)
memory usage: 33.3 KB

1. Реализовать аппроксимацию распределений данных с помощью ядерных оценок.

Распределение количественных признаков

for n in ['age', 'trestbps', 'chol', 'thalach', 'oldpeak']:
    sns.displot(x=n,
                kde=True,
                data=df)

Построены графики ядерных оценок для количесвтвенных признаков. Для качесвенных признаков нет смысла строить, так как они распределены по группам.

2. Реализовать анализ данных с помощью cdplot, dotchart, boxplot и stripchart.

Построим графики для признаков trestbps и cp, посмотрим, как распределёны значения trestbps в зависимости от значений cp.

cdplot

sns.displot(data = df,
            x='trestbps',
            hue='cp',
            multiple='fill',
            kind='kde')

<seaborn.axisgrid.FacetGrid at 0x7fb4798b0df0>

cdpolt демонстрирует разное распределение признака trestbps в зависимости от значения переменной cp.

boxplot

sns.boxplot(data=df, x='cp', y='trestbps')

<AxesSubplot:xlabel='cp', ylabel='trestbps'>

boxplot показывают неоднородность в распределении признаков.

dotchart

fig = px.scatter(df, y='cp', x='trestbps')
fig.update_yaxes(type='category', tickvals=df['cp'].tolist())
fig.show()

dotchart демонстрируют неоднородность в распределении признаков.

stripchart

plt.rcParams['figure.figsize'] = (30, 3)
sns.stripplot(data=df,
              y='cp',
              x='trestbps',
              jitter=False,
              s=15, marker="D", linewidth=1, alpha=0.1)
plt.yticks([0, 1, 2, 3])

([<matplotlib.axis.YTick at 0x7fb47907b8b0>,
  <matplotlib.axis.YTick at 0x7fb47907b130>,
  <matplotlib.axis.YTick at 0x7fb477e98940>,
  <matplotlib.axis.YTick at 0x7fb477e98970>],
 [Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, '')])

stripchart демонстрируют неоднородность в распределении признаков.

3. Проверить, являются ли наблюдения выбросами с точки зрения формальных статистических критериев Граббса и Q-теста Диксона. Визуализировать результаты.

Создадим массив с двумя выбросами и применим к нему тесты Граббса и Диксона.

chol_ = sorted(list(df.chol))
chol = chol_[140:160]
chol.append(chol_[0])
chol.append(chol_[-1])

plt.rcParams['figure.figsize'] = (20, 3)
sns.boxplot(data=chol, orient = 'h')

<AxesSubplot:>

Критерий Граббса

gtest = grubbs.test(chol)
print('Grubbs test', end=' ')
print(*gtest)
print('chol', end=' ')
print(*chol)

print('min', end=' ')
print(*grubbs.min_test_outliers(chol))
print('max', end=' ')
print(*grubbs.max_test_outliers(chol))

Grubbs test 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244
chol 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 126 564
min 126
max 564

sns.boxplot(data=gtest, orient='h')

<AxesSubplot:>

Было выявлено два значения выбросов (минимальное и максимально). Они и были удалены из выборки.

Критерий Диксона

dtest = DixonTest(dataframe=pd.DataFrame(chol), q=95, info=True, autorm=True)

Detected outlier: 
      0
21  564

print('Dixon test', end=' ')
print(*dtest[0])

Dixon test 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 126

sns.boxplot(data=dtest, orient='h')

<AxesSubplot:>

dtest = DixonTest(dataframe=pd.DataFrame(dtest[0]), q=95, info=True, autorm=True)

Detected outlier: 
      0
20  126

print('Dixon test', end=' ')
print(*dtest[0])
print('chol', end=' ')
print(*chol)

Dixon test 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244
chol 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 126 564

sns.boxplot(data=dtest, orient='h')

<AxesSubplot:>

Cначала было убрано наибольшее значение, наименьшее было убрано после второго прогона. Это произошло из-за того, что критерий Диксона за проход может убрать только одно значение.

4. Воспользоваться инструментами для заполнения пропусков в данных. Пропуски внести вручную и сравнить результаты заполнения с истинными значениями.

Сделаем 10% пропусков в данных и заполним их двумя способами.

trestbps = df[['trestbps']].copy()
trestbps.index = range(trestbps.shape[0])
missing_indexes = np.random.choice(trestbps.shape[0], int(trestbps.shape[0] * 0.1))
trestbps.loc[:, 'miss'] = trestbps.loc[:, 'trestbps']
trestbps.loc[missing_indexes, 'miss'] = np.nan

Заполним средним

mean = trestbps['miss'].fillna(trestbps['miss'].mean()).copy()
print('RMSE: ' + str(np.sqrt(np.mean((mean[missing_indexes] - trestbps.loc[missing_indexes, 'trestbps'])**2))))

RMSE: 14.606377063125159

Заполним медианой

median = trestbps['miss'].fillna(trestbps['miss'].median()).copy()
print('RMSE: ' + str(np.sqrt(np.mean((median[missing_indexes] - trestbps.loc[missing_indexes, 'trestbps'])**2))))

RMSE: 14.704874475266129

В обоих методах получили схожие значение ошибки заполнения данных.

5. Сгенерировать данные из нормального распределения с различными параметрами и провести анализ с помощью графиков эмпирических функций распределений, квантилей, метода огибающих, а также стандартных процедур проверки гипотез о нормальности (критерии Колмогорова-Смирнова, ШапироУилка, Андерсона-Дарлинга, Крамера фон Мизеса, Колмогорова-Смирнова в модификации Лиллиефорса и Шапиро-Франсия). Рассмотреть выборки малого и умеренного объемов.

seed=42
np.random.seed(seed)

#функция для генерации эталонного распределения
def generate_sample(e, sigma):
    x = np.linspace(e - 3*sigma, e + 3*sigma, 1000)
    x_ = np.linspace(-3, 3, 1000)
    y = norm.cdf(x_)
    sns.lineplot(x=x, y=y)

Сгенерируем данные из нормального распределения следующих видов: 1) N(0, 1) n = 100 2) N(0, 1) n = 1000 3) N(7, 3) n = 100 4) N(7, 3) n = 1000

N = {(0, 1, 100) : norm.rvs(loc=0, scale=1, size=100, random_state=seed),
     (0, 1, 1000) : norm.rvs(loc=0, scale=1, size=1000, random_state=seed),
     (7, 3, 100) : norm.rvs(loc=7, scale=3, size=100, random_state=seed),
     (7, 3, 1000) : norm.rvs(loc=7, scale=3, size=1000, random_state=seed)}

С помощью разных методов проверим данные на нормальное распределение.

Эпирические функции распределения

for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
    sns.ecdfplot(N[(e, sigma, n)])
    generate_sample(e, sigma)
    plt.show()

Функции распределения для разных распределений повторяют друг друга с точностью до масштаба. При увеличении количества объектов функция распрделения становиться более гладкой.

Квантили

for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
    stats.probplot(N[(e, sigma, n)], dist="norm", plot=plt)
    plt.show()

Подобно графикам функций распределения: при увеличении числа элементов выборки график становиться более гладким.

Метод огибающих

Реализация только на R.

По данному тесту выборка удовлетворяет нормальному распределению.

Процедуры проверки гипотез о нормальности

Критерий Колмогорова-Смирнова

for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
    print(f'N({e}, {sigma}), n = {n}:', kstest(N[(e, sigma, n)], stats.norm.cdf))

N(0, 1), n = 100: KstestResult(statistic=0.10357070563896065, pvalue=0.21805553378516185)
N(0, 1), n = 1000: KstestResult(statistic=0.017327787320720822, pvalue=0.9196626608357358)
N(7, 3), n = 100: KstestResult(statistic=0.9348327888393462, pvalue=5.05808847359865e-119)
N(7, 3), n = 1000: KstestResult(statistic=0.9436852219799222, pvalue=0.0)

Критерий Шапиро-Уилка

for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
    print(f'N({e}, {sigma}), n = {n}:', shapiro(N[(e, sigma, n)]))

N(0, 1), n = 100: ShapiroResult(statistic=0.9898834228515625, pvalue=0.6551706790924072)
N(0, 1), n = 1000: ShapiroResult(statistic=0.9986080527305603, pvalue=0.6264819502830505)
N(7, 3), n = 100: ShapiroResult(statistic=0.9898834228515625, pvalue=0.6551706790924072)
N(7, 3), n = 1000: ShapiroResult(statistic=0.9986080527305603, pvalue=0.6264819502830505)

Критерий Андерсона-Дарлинга

for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
    print(f'N({e}, {sigma}), n = {n}:', anderson(N[(e, sigma, n)]))

N(0, 1), n = 100: AndersonResult(statistic=0.25343395875111696, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
N(0, 1), n = 1000: AndersonResult(statistic=0.3474697767348971, critical_values=array([0.574, 0.653, 0.784, 0.914, 1.088]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
N(7, 3), n = 100: AndersonResult(statistic=0.25343395875111696, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
N(7, 3), n = 1000: AndersonResult(statistic=0.3474697767348971, critical_values=array([0.574, 0.653, 0.784, 0.914, 1.088]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))

Критерий Крамера фон Мизеса

for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
    print(f'N({e}, {sigma}), n = {n}:', cramervonmises(N[(e, sigma, n)], stats.norm.cdf))

N(0, 1), n = 100: CramerVonMisesResult(statistic=0.1549979480315905, pvalue=0.3751758581771183)
N(0, 1), n = 1000: CramerVonMisesResult(statistic=0.06313990902074608, pvalue=0.7943350162473933)
N(7, 3), n = 100: CramerVonMisesResult(statistic=31.500937472206292, pvalue=1.2769140411705848e-08)
N(7, 3), n = 1000: CramerVonMisesResult(statistic=319.7856715981252, pvalue=1.0722391730055847e-07)

Колмогорова-Смирнова в модификации Лиллиефорса

for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
    print(f'N({e}, {sigma}), n = {n}:', lilliefors(N[(e, sigma, n)]))

N(0, 1), n = 100: (0.051776473605976925, 0.7370142762533108)
N(0, 1), n = 1000: (0.021396139137066617, 0.4050632259039718)
N(7, 3), n = 100: (0.05177647360597731, 0.7370142762532998)
N(7, 3), n = 1000: (0.02139613913706656, 0.40506322590397575)

Для всех критериев значение p-value больше критического значение 0.05 у нормального рапределения N(0, 1), селедовательно принимаем гипотезу о нормальности данных. Для N(7, 3) для всех тестов кроме Колмогорова-Смирнова при n = 1000 принимаем гипотезу.

6. Продемонстрировать пример анализа данных с помощью графиков квантилей, метода огибающих, а также стандартных процедур проверки гипотез о нормальности. Рассмотреть выборки малого и умеренного объемов.

Будем рассматривать данные из колонки trestbps, так как ее график наиболее схож с графиком нормального распределения.

sns.displot(x='trestbps', kde=True, aspect=2.5, data=df)

<seaborn.axisgrid.FacetGrid at 0x7fb4763203a0>

#функция для генерации эталонного распределения
def generate_sample(e, sigma):
    x = np.linspace(e - 3*sigma, e + 3*sigma, 1000)
    x_ = np.linspace(-3, 3, 1000)
    y = norm.cdf(x_)
    sns.lineplot(x=x, y=y)

Малая выборка

Возьмём выборку из 50 элементов.

trestbps_50 = df.trestbps[0:50]

plt.figure(figsize=(18,7))

plt.subplot(1, 3, 1)
sns.ecdfplot(trestbps_50, color = 'orange')
plt.xlabel('ecdf')

plt.subplot(1, 3, 2)
generate_sample(135, 18)
sns.ecdfplot(trestbps_50, color = 'orange')
plt.xlabel('ecdf, N(135, 18)')

plt.subplot(1, 3, 3)
generate_sample(135, 18)
plt.xlabel('N(135, 18)')

plt.show()

stats.probplot(trestbps_50, dist="norm", plot=plt)
plt.show()

Метод огибающих

print(kstest(trestbps_50, stats.norm.cdf))
print(shapiro(trestbps_50))
print(anderson(trestbps_50))
print(cramervonmises(trestbps_50, stats.norm.cdf))
print('Lilliefors', lilliefors(trestbps_50))

KstestResult(statistic=1.0, pvalue=0.0)
ShapiroResult(statistic=0.9752922654151917, pvalue=0.37441587448120117)
AndersonResult(statistic=0.5554230489995859, critical_values=array([0.538, 0.613, 0.736, 0.858, 1.021]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
CramerVonMisesResult(statistic=16.666666666666664, pvalue=5.220809340400479e-10)
Lilliefors (0.10271428206231425, 0.21230434604399845)

Половина тестов указывает на то, чтобы принять гипотезу, при этом все результаты на R для малой выборки подтверждают гипотезу о нормальности распределения, значит можно говорить нормальности данных.

Умеренная выборка

Возьмём все данные признака trestbps.

trestbps = df.trestbps

plt.figure(figsize=(18,7))

plt.subplot(1, 3, 1)
sns.ecdfplot(trestbps, color = 'orange')
plt.xlabel('ecdf')

plt.subplot(1, 3, 2)
generate_sample(135, 18)
sns.ecdfplot(trestbps, color = 'orange')
plt.xlabel('ecdf, N(135, 18)')

plt.subplot(1, 3, 3)
generate_sample(135, 18)
plt.xlabel('N(135, 18)')

plt.show()

stats.probplot(trestbps, dist="norm", plot=plt)
plt.show()

Метод огибающих

print(kstest(trestbps, stats.norm.cdf))
print(shapiro(trestbps))
print(anderson(trestbps))
print(cramervonmises(trestbps, stats.norm.cdf))
print("Lilliefors", lilliefors(trestbps))

KstestResult(statistic=1.0, pvalue=0.0)
ShapiroResult(statistic=0.965917706489563, pvalue=1.4579997014152468e-06)
AndersonResult(statistic=2.5716952861943128, critical_values=array([0.569, 0.648, 0.777, 0.906, 1.078]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
CramerVonMisesResult(statistic=101.0, pvalue=0)
Lilliefors (0.10194559488884891, 0.0009999999999998899)

Все тесты указывают на то, что нужно отвергнуть гипотезу.Итого, результаты не дают нам однозначно сделать вывод о том, что данные имеют нормальное распределение.

7.

7. Продемонстрировать применение для проверки различных гипотез и различных доверительных уровней (0.9, 0.95, 0.99) следующих критериев: a. Стьюдента, включая односторонние варианты, когда проверяемая нулевая гипотеза заключается в том, что одно из сравниваемых средних значений больше (или меньше) другого. Реализовать оценку мощности критериев при заданном объеме выборки или определения объема выборки для достижения заданной мощности; b. Уилкоксона-Манна-Уитни (ранговые); c. Фишера, Левене, Бартлетта, Флигнера-Килина (проверка гипотез об однородности дисперсий).

Тест Стьюдента

Будем рассматривать данные по колонкам trestbps и thalach, так как у них более менее схожие параметры.

def student_test(alpha, data1, data2):
    # Гипотеза о равенстве
    stat, p_value = ttest_ind(data1, data2)
    if p_value > alpha:
        print('Принять гипотезу о равенстве')
    else:
        print('Отклонить гипотезу о равенстве')
    print("Мощность критерия", smp.ttest_power(np.mean(data1)/np.std(data1),
                                               nobs=len(data1), alpha = alpha, alternative='two-sided'))
    print()
    # Гипотеза о большем значении среднего первой выборке'
    stat, p_value = ttest_ind(data1, data2, alternative = 'less')
    if p_value > alpha:
        print('Принять гипотезу о большем значении среднего первой выборки')
    else:
        print('Отклонить гипотезу о большем значении среднего первой выборки')
    print("Мощность критерия", smp.ttest_power(np.mean(data1)/np.std(data1),
                                               nobs=len(data1), alpha = alpha, alternative='smaller'))
    print()
    # Гипотез о меньшем значении среднего перовой выборке
    stat, p_value = ttest_ind(data1, data2, alternative = 'greater')
    if p_value > alpha:
        print('Принять гипотезу о меньшем значении среднего первой выборки')
    else:
        print('Отклонить гипотезу о меньшем значении среднего первой выборки')
    
    print("Мощность критерия", smp.ttest_power(np.mean(data1)/np.std(data1),
                                               nobs=len(data1), alpha = alpha, alternative='larger'))

for alpha in (0.1, 0.05, 0.01):
    print('\nДоверительный интервал '+str(1 - alpha))
    student_test(alpha, df.trestbps, df.thalach)

Доверительный интервал 0.9
Отклонить гипотезу о равенстве
Мощность критерия 1.0

Отклонить гипотезу о большем значении среднего первой выборки
Мощность критерия 0.0

Принять гипотезу о меньшем значении среднего первой выборки
Мощность критерия 1.0

Доверительный интервал 0.95
Отклонить гипотезу о равенстве
Мощность критерия 1.0

Отклонить гипотезу о большем значении среднего первой выборки
Мощность критерия 0.0

Принять гипотезу о меньшем значении среднего первой выборки
Мощность критерия 1.0

Доверительный интервал 0.99
Отклонить гипотезу о равенстве
Мощность критерия 1.0

Отклонить гипотезу о большем значении среднего первой выборки
Мощность критерия 0.0

Принять гипотезу о меньшем значении среднего первой выборки
Мощность критерия 1.0

Принимаем гипотезу о том, что среднее параметра trestbps меньше среднего параметра thalach.

Тест Уилкоксона-Манна-Уитни

Рассмотрим разные выборки из признака trestbps.

def w_test(alpha, data1, data2):
    # Гипотеза о равенстве
    stat, p_value = stats.mannwhitneyu(data1, data2, alternative='two-sided')
    if p_value > alpha:
        print('Принять гипотезу о равенстве')
    else:
        print('Отклонить гипотезу о равенстве')
        
    # Гипотеза о большем значении среднего первой выборке'
    stat, p_value = stats.mannwhitneyu(data1, data2, alternative='less')
    if p_value > alpha:
        print('Принять гипотезу о большем значении среднего первой выборки')
    else:
        print('Отклонить гипотезу о большем значении среднего первой выборки')
    # Гипотез о меньшем значении среднего перовой выборке
    stat, p_value = stats.mannwhitneyu(data1, data2, alternative='greater')
    if p_value > alpha:
        print('Принять гипотезу о меньшем значении среднего первой выборки')
    else:
        print('Отклонить гипотезу о меньшем значении среднего первой выборки')

for alpha in (0.1, 0.05, 0.01):
    print('\nДоверительный интервал '+str(1 - alpha))
    w_test(alpha, df.trestbps[:30], df.trestbps[-30:])

Доверительный интервал 0.9
Принять гипотезу о равенстве
Принять гипотезу о большем значении среднего первой выборки
Принять гипотезу о меньшем значении среднего первой выборки

Доверительный интервал 0.95
Принять гипотезу о равенстве
Принять гипотезу о большем значении среднего первой выборки
Принять гипотезу о меньшем значении среднего первой выборки

Доверительный интервал 0.99
Принять гипотезу о равенстве
Принять гипотезу о большем значении среднего первой выборки
Принять гипотезу о меньшем значении среднего первой выборки

В различных выборках из признака trestbps данные схожи между собой по среднему.

Тест Фишера, Левене, Бартлетта, Флигнера-Килина

Рассмотрим признаки trestbps и thalach.

sns.displot(x='trestbps', kde=True, data=df)

<seaborn.axisgrid.FacetGrid at 0x7fb474f45700>

sns.displot(x='thalach', kde=True, data=df)

<seaborn.axisgrid.FacetGrid at 0x7fb474e6b880>

def tests(alpha, data1, data2):
    f = np.var(data1, ddof=1)/np.var(data2, ddof=1)
    x = np.array(data1)
    y = np.array(data2)
    
    nun = x.size-1
    dun = y.size-1
    p_value = 1-stats.f.cdf(f, nun, dun)
    print('\nТест Фишера:', end=' ')
    if p_value > alpha:
        print('Одно распределение')
    else:
        print('Разное рапределение')
        
    stat, p_value = stats.levene(data1, data2, center='mean')
    print('\nТест Левене:', end=' ')
    if p_value > alpha:
        print('Одно распределение')
    else:
        print('Разное рапределение')
    
    stat, p_value = stats.bartlett(data1, data2)
    print('\nТест Бартлета:', end=' ')
    if p_value > alpha:
        print('Одно распределение')
    else:
        print('Разное рапределение')
               
    stat, p_value = stats.fligner(data1, data2, center='mean')
    print('\nТест Флигнера-Килина:', end=' ')
    if p_value > alpha:
        print('Одно распределение')
    else:
        print('Разное рапределение')

for alpha in (0.1, 0.05, 0.01):
    print('\nДоверительный интервал '+str(1 - alpha))
    tests(alpha, df.trestbps[:30], df.thalach[:30])

Доверительный интервал 0.9

Тест Фишера: Одно распределение

Тест Левене: Одно распределение

Тест Бартлета: Одно распределение

Тест Флигнера-Килина: Одно распределение

Доверительный интервал 0.95

Тест Фишера: Одно распределение

Тест Левене: Одно распределение

Тест Бартлета: Одно распределение

Тест Флигнера-Килина: Одно распределение

Доверительный интервал 0.99

Тест Фишера: Одно распределение

Тест Левене: Одно распределение

Тест Бартлета: Одно распределение

Тест Флигнера-Килина: Одно распределение

Признаки trestbps и thalach имеют одно и то же распределение.

8. Исследовать корреляционные взаимосвязи в данных с помощью коэффициентов корреляции Пирсона, Спирмена и Кендалла.

Построим корреляционные матрицы.

for c in ['pearson', 'spearman', 'kendall']:
    plt.figure(figsize = (24, 12))
    sns.heatmap(df.corr(method=c),
                annot=True,
                fmt=".2f")
    plt.title(c)
    plt.show()

Признаки очень слабо коррелируют между собой. Самая высокая по модулю корреляция между признаками slope и oldpeak.

9. Продемонстрировать использование методов хи-квадрат, точного теста Фишера, теста МакНемара, Кохрана-Мантеля-Хензеля.

Метод хи-квадрат

Проверим независимость колонок sex (пол) и target (наличие болезни).

con = pd.crosstab(df.sex, df.target)

ret = chi2_contingency(con)[1]
print(round(ret, 5))

0.0

Отвергаем нулевую гипотезу о независимости пола и наличия болезни.

Точный тест Фишера

Проверим независимость колонок sex (пол) и target (наличие болезни).

ret = stats.fisher_exact(con)[1]
print(round(ret, 5))

0.0

Отвергаем нулевую гипотезу о независимости пола и наличия болезни.

Тест МакНемара

Проверим независимость колонок sex (пол) и target (наличие болезни).

print(mcnemar(np.array(con), exact=False)) #используем распределение хи квадрат

pvalue      0.00264477649712634
statistic   9.03763440860215

Отвергаем нулевую гипотезу о независимости пола и наличия болезни.

Тест Кохрана-Мантеля-Хензеля

Проверим независимость колонок sex (пол), target (наличие болезни) и exang (стенокардии, вызванной физической нагрузкой).

data = df[['sex', 'exang', 'target']]

CMH(data, 'sex', 'target', stratifier = 'exang')

Cochran-Mantel-Haenszel Chi² test

sex x target, stratified by exang

Cochran-Mantel-Haenszel M^2 = 18.04335, dof = 1, p-value = 0.0

Следовательно отвергаем нулевую гипотезу о независимости пола, наличии болезни и стенокардии, вызванной физической нагрузкой.

10. Проверить наличие мультиколлинеарности в данных с помощью корреляционной матрицы и фактора инфляции дисперсии.

Корреляционная матрица

Построим корреляционную матртицу

df.corr()

	age	sex	cp	trestbps	chol	fbs	restecg	thalach	exang	oldpeak	slope	ca	thal	target
age	1.000000	-0.098447	-0.068653	0.279351	0.213678	0.121308	-0.116211	-0.398522	0.096801	0.210013	-0.168814	0.276326	0.068001	-0.225439
sex	-0.098447	1.000000	-0.049353	-0.056769	-0.197912	0.045032	-0.058196	-0.044020	0.141664	0.096093	-0.030711	0.118261	0.210041	-0.280937
cp	-0.068653	-0.049353	1.000000	0.047608	-0.076904	0.094444	0.044421	0.295762	-0.394280	-0.149230	0.119717	-0.181053	-0.161736	0.433798
trestbps	0.279351	-0.056769	0.047608	1.000000	0.123174	0.177531	-0.114103	-0.046698	0.067616	0.193216	-0.121475	0.101389	0.062210	-0.144931
chol	0.213678	-0.197912	-0.076904	0.123174	1.000000	0.013294	-0.151040	-0.009940	0.067023	0.053952	-0.004038	0.070511	0.098803	-0.085239
fbs	0.121308	0.045032	0.094444	0.177531	0.013294	1.000000	-0.084189	-0.008567	0.025665	0.005747	-0.059894	0.137979	-0.032019	-0.028046
restecg	-0.116211	-0.058196	0.044421	-0.114103	-0.151040	-0.084189	1.000000	0.044123	-0.070733	-0.058770	0.093045	-0.072042	-0.011981	0.137230
thalach	-0.398522	-0.044020	0.295762	-0.046698	-0.009940	-0.008567	0.044123	1.000000	-0.378812	-0.344187	0.386784	-0.213177	-0.096439	0.421741
exang	0.096801	0.141664	-0.394280	0.067616	0.067023	0.025665	-0.070733	-0.378812	1.000000	0.288223	-0.257748	0.115739	0.206754	-0.436757
oldpeak	0.210013	0.096093	-0.149230	0.193216	0.053952	0.005747	-0.058770	-0.344187	0.288223	1.000000	-0.577537	0.222682	0.210244	-0.430696
slope	-0.168814	-0.030711	0.119717	-0.121475	-0.004038	-0.059894	0.093045	0.386784	-0.257748	-0.577537	1.000000	-0.080155	-0.104764	0.345877
ca	0.276326	0.118261	-0.181053	0.101389	0.070511	0.137979	-0.072042	-0.213177	0.115739	0.222682	-0.080155	1.000000	0.151832	-0.391724
thal	0.068001	0.210041	-0.161736	0.062210	0.098803	-0.032019	-0.011981	-0.096439	0.206754	0.210244	-0.104764	0.151832	1.000000	-0.344029
target	-0.225439	-0.280937	0.433798	-0.144931	-0.085239	-0.028046	0.137230	0.421741	-0.436757	-0.430696	0.345877	-0.391724	-0.344029	1.000000

Признаки очень слабо коррелируют между собой. Самая высокая по модулю корреляция между признаками slope и oldpeak.

Фактор инфляции дисперсии

Применим данный метод ко всем признакам.

y, x = dmatrices('target ~ age + sex + cp + trestbps + chol + fbs + restecg + thalach + exang + oldpeak + slope + ca + thal', data=df, return_type='dataframe')
vif = pd.DataFrame()
vif['VIF'] = [variance_inflation_factor(x.values, i) for i in range(x.shape[1])]
vif['variable'] = x.columns
vif

	VIF	variable
0	207.256646	Intercept
1	1.443474	age
2	1.161866	sex
3	1.284456	cp
4	1.170591	trestbps
5	1.150174	chol
6	1.087379	fbs
7	1.060998	restecg
8	1.613726	thalach
9	1.402001	exang
10	1.705857	oldpeak
11	1.642595	slope
12	1.202570	ca
13	1.147279	thal

Получились значения близкие к единице, что указывает на очень слабую коллинеарность признаков.

11. Исследовать зависимости в данных с помощью дисперсионного анализа.

Применим данный метод ко всем признакам в паре с target.

for pr in ['age', 'sex', 'cp', 'trestbps', 'chol', 'fbs', 'restecg', 'thalach', 'exang', 'oldpeak', 'slope', 'ca', 'thal']:
    aov = stats.f_oneway(df[pr], df.target)[1]
    print(pr+' ~ target p-value =', aov)

age ~ target p-value = 0.0
sex ~ target p-value = 0.0004397191734875478
cp ~ target p-value = 2.845231775397828e-10
trestbps ~ target p-value = 0.0
chol ~ target p-value = 0.0
fbs ~ target p-value = 8.910856675010132e-27
restecg ~ target p-value = 0.6920275427852707
thalach ~ target p-value = 0.0
exang ~ target p-value = 4.686916079913056e-08
oldpeak ~ target p-value = 2.2270855709918155e-11
slope ~ target p-value = 2.970631656050962e-62
ca ~ target p-value = 0.004846478356075194
thal ~ target p-value = 4.408516795472436e-167

Отвергаем нулевую гипотезу об отсутствии влияния признака на наличие болезни для age, sex, cp, trestbps, thalach, exang, oldpeak, ca, thal, chol, fbs, slope; Принимаем нулевую гипотезу для restecg.Для R получены другие результаты для chol, fbs, slope.

12. Подогнать регрессионные модели (в том числе, нелинейные) к данным, а также оценить качество подобной аппроксимации.

Попробуем предугадать значение trestbps (давление в состоянии покоя) от age (возраста).

def regr_model(X, y):
    
    X = np.array(X)[:, np.newaxis]
    
    regr = LinearRegression()
    cub = PolynomialFeatures(degree=3)
    
    X_fit = np.arange(X.min(), X.max(), 1)[:, np.newaxis]
    X_cub = cub.fit_transform(X)
    
    regr = regr.fit(X, y)
    y_lin_fit = regr.predict(X_fit)
    
    regr = regr.fit(X_cub, y)
    y_cub_fit = regr.predict(cub.fit_transform(X_fit))

    plt.scatter(X, y)

    plt.plot(X_fit, y_lin_fit, 
             label='linear',
             lw=2)

    plt.plot(X_fit, y_cub_fit, 
             label='cubic',
             lw=2)

    plt.legend()

plt.rcParams["figure.figsize"] = (6, 4)
regr_model(df.age, df.trestbps)

Построили две модели, получилось не очень хорошо (судя по точкам на графике и не могло получится хорошо). Наверно, влияет, что в выборке люди как с болезнью, так и без.

Вывод

-Узнал о 12 статистических методах -Научился их применять с помощью языков Python и R -На R удобнее реализованы данные методы -Некоторые методы дают разные результаты в разных языках на одних и тех же данных -Интересно и полезно провёл время